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Time-dependent quantum mechanics provides an intuitive picture of particle propaga- 
tion in external fields. Semiclassical methods link the classical trajectories of particles with 
their quantum mechanical propagation. Many analytical results and a variety of numerical 
methods have been developed to solve the time-dependent Schrodinger equation. The time- 
dependent methods work for nearly arbitrarily shaped potentials, including sources and sinks 
via complex- valued potentials. Many quantities are measured at fixed energy, which is seem- 
ingly not well suited for a time-dependent formulation. Very few methods exist to obtain the 
energy-dependent Green function for complicated potentials without resorting to ensemble av- 
erages or using certain lead-in arrangements. Here, we demonstrate in detail a time-dependent 
approach, which can accurately and effectively construct the energy-dependent Green func- 
tion for very general potentials. The applications of the method are numerous, including 
I chemical, mesoscopic, and atomic physics. 

1 Introduction 

In principle, the energy Green function G(r,r';_E) contains all one ever needs to know about a 
stationary system. For example, in atomic systems, the Green function gives the photoionization 
and photodetachment rate in the presence of electric and magnetic fields [TSJ [3]- The theory 
of atom lasers emitted from Bose-Einstein condensates relies on the Green function [16 . The 

• *h Green function is of course quite nontrivial to obtain, even for single particle problems, if disorder, 

/\ magnetic fields, etc. are present. 

In mesoscopic systems (which are characterized by large free path lengths) the Green function 
connects the transport of electrons through the system with the conductivity, and makes possible 
the determination of the microscopic electron flow. With the advent of imaging techniques which 
spatially resolve the electron flow on the nano scale, the realistic modelling of a device, including 
charged gates, disorder potentials and screening, has become a necessity |20j . 

Traditionally, Green functions for mesoscopic systems demanded either a diagrammatic ap- 
proach based on perturbation theory (chap. 5 of [?]) or to connect leads to the scattering region 
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to obtain the Green function of the full system consisting of leads and scattering region by re- 
cursion (chap. 3 of 8 ). The diagrammatic approach has the serious drawback that it is limited 
to the ensemble average of Green functions for many realizations of the potential. The ensemble 
averaged current through a device is in general very different from the flow through a typical 
realization of the device. The Green function for a single realization of the device, instead of an 
ensemble, has become central to many experiments. The recursive Green function technique has 
serious limitations due to the attachment of waveguide-type leads to the system which have to be 
arranged colinearly with the scattering region (the last condition can be relaxed by employing a 
smart reordering scheme based on graph theory 21J). 

The approach presented in this paper permits us to obtain the Green function for a single real- 
ization with high accuracy and short computational time - without the introduction of leads. Our 
approach works for very general electron-sources attached to a device, i.e. the injection of current 
to the two-dimensional electronic subsystem from a (3D) reservoir through a tunneling contact 
(somewhat similar to the current transport from the tip of a scanning tunneling microscope [B]). 
Also inner contacts (which are used i.e. in Mach-Zehnder interferometers in mesoscopic systems) 
are possible to incorporate in the time-dependent method, as well as spin-orbit interactions. 

The Green function contains the discrete cigcncncrgics E n and the continuous ones E v , and 
all eigcnfunctions ip n ,u of the Hamiltonian in a compact form: 

G(r y ;£0=E M^) + j^^mMi. (d 

Surprisingly, only very few energy-dependent Green functions are known analytically, in contrast 
to a vast amount of available time-dependent propagators [TU] . The direct calculation of the Green 
function for a given potential is generally not possible and perturbation expansions diverge. Thus 
the spectral properties of many systems arc unknown, even though they are very relevant for 
predicting and analyzing experimental results. 

One way of obtaining an energy Green function is via Fourier transform from the time domain[18 
1191 1131 [T2"llll) . This can be done numerically, as it is always a one dimensional integral no matter 
how many degrees of freedom are present. However this method of obtaining the Green function 
has not been fully exploited or explored. In this paper, we describe an approach for obtaining 
energy-dependent Green functions via the time-dependent Schrodinger equation. In Sect. [2] we 
establish the connection between the time-dependent Schrodinger equation, the autocorrelation 
function of a wavepacket, and the energy-dependent Green function. Analytical and numerical 
examples of the method are discussed in Sects. [3] and |4j We provide a quantitative analysis of the 
accuracy of the method in Sect. [5j We apply our method to electron propagation through a nano 
device in Sect. [6] In Sect. [7] we discuss thermal wavepackets and finally give some conclusions. 



2 The time- dependent Schrodinger equation. 

In the time-dependent formulation of quantum mechanics, the use of the propagator K permits us 
to transform the partial differential time-dependent Schrodinger equation for an initial quantum 
state ip(r',ti) into an integral equation: 

fl>(T,t f ) = f dr' K{r,t f \r',ti)i>{r',ti) (2) 

For a time-independent Hamiltonian, we can set t = tt—ti. The energy-dependent Green function 
is tied to the time-dependent propagator K via a Laplace transform: 

G(r, r' ;£) = - lira / dt S E+i ^ t / h K(r, t\r' , 0), (3) 
ih v ^o+ J 

where rj is a small positive number, which selects the outgoing wave boundary condition at r and 
thus the retarded Green function. However, the Laplace transform is seldomly used to obtain the 
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Green function analytically or numerically (even if the propagator is known), since the propagator 
is a highly oscillatory function with additional singularities. 

For Hamiltonians which are polynomials of positions and momenta up to quadratic order, the 
time-dependent propagator assumes the form 

K(r,t\r',t') = a(t)exp[iS(r,t\r',t')/K], (4) 

where S denotes the classical action, which is itself a maximally quadratic polynomial in positions. 
Therefore it is possible to integrate Eq. ^ analytically for Gaussian wavepackets in D dimensions 
since all integrals involve Gaussian functions. For non-quadratic Hamiltonians, the analytic form 
of the propagator is in general not known. However, various methods have been developed to solve 
the time-dependent Schrodinger equation numerically (see Sec. [4]). 

2.1 From wavepackets to the Green function 

The energy-dependent Green function satisfies the inhomogeneous Schrodinger equation 

[H-E]G{r,r';E)=S(r-r') (5) 

with a ^-distribution "source" -term. If we specify as source-term a spatially extended wavepacket 
ip(r), we obtain the overlap of the Green function with the state 

(i>\G{E)\i>) = J dr J dr'i[>(r,0)*G(r,r';E)>ip(r',0) 

= ^J™ dt eiEt/h / dr / dr '^ (r ' °)*^( r ' *l r '' °)^( r '' °) 
= ijf dte lEt ' h J drV>(r,0)*V>(r,i) 

dte iEt/h C(t), (6) 



m 

where C(t) denotes the autocorrelation function 

C(t) = J dr^(r,Q)*V(r,i). (7) 

The Green function for r' ^ r is obtained by replacing (ip\ with another state (tp r i | centered around 
r' while keeping \ip). In numerical calculations r' is a grid point of the discretized representation 



of the potential and of the wave functions (see eq. (35 1) 



2.2 Calculation of the local density of states (LDOS) 

The local density of states (LDOS) at point r is the imaginary part of the Green function 

n(r;E) = SG(r,r;^. (8) 

7T 

The definition of the LDOS requires to use a <5-distributed initial state, which cannot be represented 
by a finite width wavepacket. This is an important technical point, requiring some care to handle 
correctly. Essentially, a point source (as in the coordinate space Green function) implies very 
high (in fact infinite) energy components, which are inconvenient and physically unimportant in 
a practical sense. In Sect. [5] we discuss how to recover the Green function from a finite width 
wavepacket propagation, to very high accuracy. 
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3 Analytical Examples: Free particle and motion in a uni- 
form electric field 



In the following, we take as initial wavepacket in two dimensions a Gaussian of the form 



V>(r) 



: exp 



r 

2d 2 



(9) 



The probability density of this wavepacket is normalized to unity. We obtain the LDOS in correct 
units, if we multiply the correlation function by 




1 



(10) 



where the first contribution is the switch of the normalization from unity of the modulus squared 
J |-0(r)| 2 d 2 r = 1 to the non-squared function J -0(r)d 2 r = 1, (this in order to correspond to the 
delta function normalization of a coordinate state |r)) and the second factor originates from getting 
the imaginary part by using Laplace transform . Thus we obtain for the wave-packet-LDOS for a 
wavepacket centered around r 



,(«) 



(r;E)=- 



1 

47r 2 a 2 



dt e iEt/n C(t) 



(11) 



which depends on parameters of the wavepacket (like the width a) . For numerical reasons we do 
not want to be propagating wavepackets with very high energy components if we can avoid it. In 
order to discover how to compensate for the finite width of the wavepacket, we calculate specific 
autocorrelation functions analytically for systems with continuous spectra. 



3.1 The free propagator in 2D 

For the Hamiltonian 



the propagator reads 



fffree = (12) 

2m v ' 



The four Gaussian integrations to obtain the autocorrelation function starting with the wavepacket 
in Eq. ^ are straightforward and yield 

c -«> - dSW < 14) 

The Laplace transform is given by 



4n 2 a 2 Air 2 a 2 iTi 



r°° p --ia mJi/h m 

/ dtC hcc (ty Et/h - —2 [27ri + T(0,-2a 2 mE/h 2 )] , (15) 

Jo 2tt z Ti 



where T(a,x) denotes the incomplete Gamma function. Using QT(0, — |x|) = — 7r, we obtain the 
LDOS for the wavepacket centered around the origin o 



n£l(o ;J B)=9 



47r 2 a 2 



™ e -^ 2 mE/H\ (16) 

2tt% 2 
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The exponential decay of the wavepacket LDOS reflects the probability density of the kinetic 
energy of the initial wavepacket p(-Ekm)i which suppresses higher energies exponentially. In the 

momentum representation, with k = + fc 2 , , the wavepacket becomes 

Mk,0) = JL / dr 2 e irk ^(r,0) = 2a^c-^ k \ (17) 
2nn J 

and the energy probability distribution p(-Ekin) = 27rm|-0(fc = t/2rnE^/'h, 0)| /fi 2 is properly 
normalized 

/ p(E kin )dE kin = 2n dk\i>{k)\ =1. (18) 
Jo Jo 

For a = we obtain the LDOS for a ^-distribution. 



3.2 The linear potential in 2D 

For the linear potential with Hamiltonian and propagator 

P 2 x+pI „ , , , m Am|r-r'| 2 \Ft(x + x') iF 2 t 3 \ 

H ™ = ^T- Fx > K ^ V> 0) = 2^i exp y 2ht + "V^ - 24^ ) 

we obtain the following autocorrelation function 



(19) 



_ , . 2a 2 m ( a 2 F 2 t 2 F 2 t 3 \ , , 

= 2^ 6XP \~^r - 1 2mn)- (20) 



The LDOS for a Gaussian wavepacket of width a centered around the origin o becomes 



T) (a) (o- E) - m 32a 6 /3-8i;a 2 /3 
"-field 1 -"' a ) — *2 e 
27T?l 



with 

1/3 



i-Aii (2 2 / 3 (4a 4 - 2£/3) 



(21) 



TO \ /J Z" 2 

'=(^sj , a = /3Fa, M i( z ) = J Q dxAi(x) (22) 



For a — > we recover the LDOS for a ^-distribution 



nfioid(o; £/) = — — 2 
27rn, 



i-Ai! (2 2 / 3 (-2£/3)) 



(23) 



4 Numerical method for general potentials 



The time evolution operator U for a time-independent Hamiltonian H with kinetic energy T and 
potential energy V can be divided into N steps: 



U(t, t') = e-Wit-t^/H _ c -i(T+V)NAt/h _ f-i(T+V)At/h\ ^ 

where At = (t — t')/N. We separate the kinetic and potential energy terms using the Baker- 
Campbell-Hausdorff series 

e -iAt/fi (V/2+T+V/2) _ e -iAt/fi V/2 e -iAt/h Tg-iAt/fi V/2 e O(At 3 ) ^5) 

The symmetric distribution of V/2 cancels all terms of order At 2 and makes the approximation 
accurate up to order At 3 [7J. The application of the kinetic energy operator e ~ lAt / h T to a wave 
function reduces to a multiplication if the wave function is given in momentum space, whereas the 
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Figure 1: Local density of states in a homogeneous electric field. Left panel: LDOS for a 6- 
distribution (a = 0), and for a — 10 nm. Also shown is Af rcc (_E). Right panel: Absolute error in 
units of the free-particle LDOS. Parameters: F — 50 keV/m, effective mass m = 0.067 m e . 



potential energy operator is a multiplication in position space. We denote the Fourier transforms 
between position and momentum space by T: 

^ (r) " pdjw / d " e " PV(t) ' ^ V(P) = / dp e ""* (p) - (26) 

Combining N consecutive steps and joining the half steps within the loop yields the solution of 
the time-dependent Schrodinger equation 

iP(r,t' + NAt) ~ e ^/nv/2 

Numerically the wave function is represented on a grid and the value of ip on the grid points 
is used to perform a Discrete Fast Fourier Transform (DFFT). Since we apply the numerical 
wavepacket method also to systems where the wavepacket will eventually leave the initial region 
of the numerical grid-representation, we have to either move the Fourier grid along with the center 
of the wavepacket or to remove the wavepacket at the borders of the region of interest. Otherwise 
the DFFT algorithm will let the wavepacket reappear at the other edge of the grid. The removal 
of the wavepacket in certain regions is achieved best by using a complex valued potential. By 
tracking the probability density, the absorbing potential region can actually simulate i.e. a contact 
(acting as an electron sink) and measure the electronic current taken out of the system in that 
region. 

5 Extracting the (5-distribution LDOS from the autocorre- 
lation function 



T7-l e -iAt/R T 



-iAt/ft V 



iAt/ft V/2 



(27) 



With the examples above in mind we can see how to compensate for the energy cutoff implied by 
a finite width wavepacket. For a wavepacket with finite width a, the energy-dependent spectrum 
is exponentially damped for energies differing greatly from E w 2} " q2 . Choosing a smaller width a 
is theoretically possible, but this would impose serious limitations for the numerical propagation 
of wavepackets, demanding a finer grid of points in coordinate and smaller time steps. The 
wavepacket LDOS has to be corrected for comparison with the ^-distribution LDOS. The correction 
stems from the initial probability distribution of the kinetic energy. For the free-particle, the 



correction factor is given by Eq. (18 1 



(a) / 

A ircc (E) = nfrec r; = exp (-2a 2 mE/h 2 ) . (28) 
n froc (r; £) 
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Figure 2: Two-dimensional potential landscape 
for magnetic focusing. The T-shaped potential of 
two quantum point contacts forms two constric- 
tions, while the background charges produce an 
irregular potential landscape. 




Energy [meV] 

Figure 3: Probability distribution of the kinetic 
energy in the initial wavepacket. Parameters: 
a = 20 nm, k = y/2mE /h, E = 9 meV, 
effective mass m = 0.067m e . 



The correction factor is dependent on the Hamiltonian under consideration, but Af rco (£') provides 
a lowest order estimate. In general, the correction factor is chosen to reproduce the kinetic energy- 
distribution of the initial wavepacket. 

The analytic results of the previous section allow to analyze the error due to the use of the 
field-free correction factor for the linear electric field (see Fig. [T]). The LDOS for a = and 
a = 10 nm are shown in the left panel. The right panel shows the absolute error of the corrected 
finite width calculation compared with the a = result. For the chosen values, the relative error 
is approximately 1/100. 

6 Application: Magnetic focusing 

The inclusion of magnetic fields necessitates splitting up the kinetic energy into two parts and 
to using a mixed momentum and position space representation. The use of a complex potential 
V a bs(r) is of special importance in the presence of magnetic fields, where periodic boundary con- 
ditions often fail to model the physical system under consideration. Magnetic focusing [2] is such 
a case; here, periodic boundary conditions result in spurious edge currents around the chosen unit 
cell, which are absent in the real device due to additional contacts. The magnetic focusing setup 
consists of two separated constrictions (quantum point contacts) in a T-shaped two-dimensional 
device. The potential is shown in Fig. [2} The irregular potential landscape above the T-shaped po- 
tential produced by the gates is due to the presence of only partially screened background charges. 
It is important to model the effective potential landscape of the device realistically. Our model 
includes the effects of screening of gate and donor atom potentials and thus includes interactions 
on a mean field level [§1 El |S] . Experimentally, Ohmic contacts before the right and after the 
left constriction are connected to different voltages, which results in a current flowing through 
the system. The whole device is placed in a homogeneous magnetic field perpendicular to the 
two-dimensional plane. 

The vector potential for a homogeneous magnetic field B = (0, 0, B) in the symmetric gauge 
becomes A = (— y, x, 0)B/2. The Hamiltonian splits into three parts 

H= p L -^ L p x y+ + W£ityc+iw2 (x 2 + y 2 ) + V(x,y), u> L =~, (29) 
2m 2m 2 2m 

■> v ' v „ ' v v ' 

T Px,y T v y ,x VeS{x,y) 

where the mixed momentum-position representation for the kinetic energy is possible since \p x , y] = 
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Figure 4: Current flow through the device at E = 9 meV. The wavepacket originates from r' = 
(0, 0). A magnetic field of B — 0.8 T is applied perpendicular to the x — y-plane. The dark regions 
denote an enhanced probability |G(r nm , r'; E)\ , eq. (35 1. The arrows denote the local current 
density, eq. ( 36 1 . The right panel is a close-up of the current flow around an obstacle placed at 
r = (—350, 140) nm. Parameters: Effective mass m = 0.067 m e . 



[p y ,x] — 0. The new propagation algorithm becomes 



V>(r, t' + NAt) 



,-iAt/ft V B ff/2 



p -iAt/h V abs T -l -iAt/h T PViX -p TT-lp-iAt/fi T px , v <p „-iAt/S V e f 



N 



xe iAt / RV ° f >/ 2 iP(r,t'), 



(30) 



where T x ,J- y denote partial Fourier transforms with respect to only one-dimension. 

In the following, we apply the wavepacket method in order to obtain a detailed picture of the 
current flow through the device. At t — we start a wavepacket which represents an expanding 
"smoke ring" centered around the origin : 



ip(k x ,k y ,t = 0) oc 



2tt 



exp 



^— a 2 /2[(k x — k Q cos 8) 2 + (k y — k sin By 



(31) 



In radial coordinates r = \J x 1 + y 2 , k — \ k 2 + k 2 the properly normalized wavefunctions read 



,J (fc r), ^(fc,t = Q) 



-a 2 k 2 /2-l/4a 2 k„ 



a^/ir I (a 2 fc 2 /2) " " ^ h(a 2 k 2 /2) 

where Jo and Io denote Bessel functions [1]. The kinetic energy probability density 

2jnn $(k = V2mE/Ti,t = 0) 



I (a 2 kk ), (32) 



(33) 

is shown in Fig. |3] The choice of a and ko is such, that the energy range under consideration is well 
covered. The wavepacket propagates via the DFFT method in N = 12000 steps of AT = 10~ 15 s. 
The grid consists of 300 x 300 grid cells of width 2.67 nm. The propagation takes about five 
minutes on a personal computer (2 GHz CPU) . Along the four edges of the grid we have placed a 
complex- valued potential of the form 



V a bs{n) = i Vo.abs cosh" 



{t% '"edge) 

J 2 



(34) 



with Vb,abs = 25 meV, d — 25 nm. During the propagation we track for each grid point r n> . 
several quantities, including 



G(r„, m ,r';£) = 



JV 



1 



A(E) 



(35) 
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Figure 5: A thermal wavepacket emerging from the quantum point contact (QPC) encounters two 
impurities which are approximately at the same radial distance from the QPC; the amplitude from 
these can interfere at the QPC (right); the amplitude scattered from impurities encountered later 
cannot interfere with the amplitude from the first two, since the arrival of that amplitude will be 
too late. 

for a set of energies E. If we want to obtain the Green function for a reasonable small region (i.e. 
around the two constrictions), we can store ip(r n<m , jAt) and obtain the energy-dependent Green 
function via the Laplace transform for any energy we pick. Thus a single run of the wavepacket 
yields G(r„ m ,r'; E) for many energies and facilitates i.e. thermal averaging considerably. 
In Fig. [4] we display the microscopic current density originating from r' 

j(r„, m , r'; E) = -3 [G(r„. m , r'; £)*VG(r n , ra , r'; E)] + -A\G(r n>m , r';E)\ 2 (36) 
m m 

throughout the whole device for E = 9 meV. The obstacle around r = (—350, 140) nm deflects the 
current and lowers the transmission from the right constriction to the left one. Experimentally, 
it is possible to use a moveable charged tip as obstacle which is moved over the sample. High- 
resolution spatial information (10 nm) of the electronic current is obtained by recording the change 
in the transmitted current between the fixed constrictions as a function of tip position [2] . The 
resulting pictures reveal interference structures and flux bundles, along which the electronic current 
is flowing through the device. 

7 Wavepackets mimicking a thermal energy distribution 

In order to obtain the energy-dependent Green function, we have discussed ways to largely remove 
the effect of wavepacket parameters, like the halfwidth a of the initial wavepacket. However, it 
is also possible to retain the wavepacket properties and use them to model a certain occupation 
of initial energies. In the case of an adiabatic quantum point contact (QPC) with one open 
transverse mode, a thermal "smokcring" wavepacket with an energy distribution matching the 
derivative of the Fermi-Dirac distribution —dfFD(E', Ep, T)/dE' replaces thermal averaging [13] . 
It is interesting that higher temperatures are represented by more confined wavepackets in position 
space, which reduces interference effects due to the shorter wave train. When considering different 
modes of backscattering which reduce the net flux through the QPC, the question arises as to 
whether they can interfere. The answer is that they interfere only if the backscattered flux from 
the thermal wavepacket arrives at the QPC at the same time (see Fig. [5]). This is a form of "white 
light interferometry" . 
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The thermal wavepacket approach permits a powerful intuitive picture and a simple estimation 
tool for understanding interference effects and electron choreography in small devices. The basic 
scenario is as follows: Launched from just to the left of the QPC, a spatially narrow wavepacket 
emerges to the right, fanning out as part of an annulus, then suffering small angle scattering, 
splitting up into pieces due to collisions with large and small objects, walls, etc. (Fig. pp. The 
pieces will have the same width in their direction of travel as the original wavepacket (unless the 
scattering is resonant and therefore time delaying). As the pieces arrive at the original QPC (or 
another terminal) the rule is simple: separate pieces of the wavepacket must arrive at the same 
time if they are to interfere. (If the pieces exit to a lead with several modes open, they interfere 
mode by mode). If an object which the wavepacket encounters is moveable, we may follow the 
oscillation (in conductance) as the interference of the objects' phase changes from constructive to 
destructive and back to constructive as it is moved in certain directions. 

8 Conclusions 

We have obtained the energy-dependent Green function for general potentials via time-dependent 
wavepacket methods. The analytical solved cases allow to quantify and to minimize the error of 
the method due to the influence of wavepacket parameters. The wavepacket approach is especially 
useful for setups where periodic boundary conditions lead to non-physical results. The modelling 
of sources and sinks by complex potentials permits us to use very realistic potentials for devices 
and contacts and is at the same time efficient and fast. The time-dependent approach allows to 
postselect the energy (one or several) for which the energy-dependent Green function is needed. 
In mesoscopic physics, the wavepacket Green function describes transport in realistic nano device 
potentials and leads to new insights for the design and operation of future devices. A straightfor- 
ward extension of the wavepacket methods allows one to include spin-orbit interactions by using 
two wavepackets corresponding to the two spin-components. 
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